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A robust rocket engine combustor design and development process must 
include tools which can accurately predict the multi-dimensional thermal 
environments imposed on solid surfaces by the hot combustion products. Currently, 
empirical methods used in the design process are typically one dimensional and do 
not adequately account for the heat flux rise rate in the near-injector region of the 
chamber. Computational Fluid Dynamics holds promise to meet the design tool 
requirement, but requires accuracy quantification, or validation, before it can be 
confidently applied in the design process. 

This effort presents the beginning of such a validation process for the Loci- 
CHEM CFD code. The model problem examined here is a gaseous oxygen 
(G02)/gaseous hydrogen (GH2) shear coaxial single element injector operating at a 
chamber pressure of 5.42 MPa. The G02/GH2 propellant combination in this 
geometry represents one the simplest rocket model problems and is thus 
foundational to subsequent validation efforts for more complex injectors. 

Multiple steady state solutions have been produced with Loci-CHEM 
employing different hybrid grids and two-equation turbulence models. Iterative 
convergence for each solution is demonstrated via mass conservation, flow variable 
monitoring at discrete flow field locations as a function of solution iteration and 
overall residual performance. A baseline hybrid grid was used and then locally 
refined to demonstrate grid convergence. Solutions were also obtained with three 
variations of the k-omega turbulence model. 

Background 

Design of reliable and durable combustion devices will be a key element in achieving 
the NASA Vision for Space Exploration. This analysis effort in the injector/chamber 
thermal compatibility area focuses on validating Computational Fluid Dynamics (CFD) 
tools for injector and chamber design by demonstrating accuracy levels for prediction of 
injector-generated thermal environments in the combustion chamber. 

Almost 50 years of experience has proven that injector/combustor failures are 
typically the result of locally severe environments and, thus, multi-dimensional. As an 
example, note the severe injector-induced blanching and cracking in a main combustion 
chamber wall as shown in Figure 1 . Given what are sure to be aggressive schedules and 
tight budgets for Space Exploration engine development, new injector design tools are 
needed to help reduce both the development time and cost of the combustion devices. 
Combustor thermal environments are a direct consequence of injector design and 
operation. Current injector design tools are mostly one-dimensional and empirical. The 
design tool shortcomings force either component over-design or an extensive, costly test 



program to reveal and fix serious design issues. Accurate, efficient tools are essential for 
the design of robust combustion devices hardware. 



Figure 1. Blanching and cracking — evidence of multi-dimensional environments in a 

combustion chamber. 

CFD holds great, but to date, largely unfulfilled promise for predicting the 
complex, three-dimensional environments during the design cycle time frame. Currently, 
injector designers continue to use their legacy tools despite the well-known issues 
associated with their use. They correctly point out that the CFD solutions are too time- 
consuming and are of questionable accuracy. The realization of CFD’s promise as an 
injector design tools will require the demonstrated improvement in the three design tool 
requirements: simulation fidelity in terms of physics, boundary conditions and geometry; 
simulation robustness in terms of ability to produce solutions in the design cycle time 
frame; and simulation accuracy in terms of comparison to relevant data. 

The combustion CFD technology effort at NASA/Marshall Space Flight Center is 
guided by a Combustion Devices CFD Simulation Capability Roadmap 1 . The Roadmap 
objective is to enable the use of CFD as a tool for the Simulation of Preburners, Ducting, 
Thrust Chamber Assemblies and Supporting Infrastructure in terms of Performance, Life, 
and Stability so as to affect the design process in a timely fashion. 

If CFD is to be used as an injector design tool, code developers & code users must 
address the following key issue. How should confidence (i.e. demonstrated accuracy 
capability) in simulations and modeling for design be critically addressed, and where 
necessary, improved? Verification and Validation of computational solutions are the 
primary means to quantify and build this confidence. 

The definitions 2 of Verification and Validation as applied to CFD analysis bear 
repeating. Verification is the process of determining that a model implementation 
accurately represents the developer’s conceptual description of the model and the 
solution to the model. Significant verification efforts 3 for the CFD tool applied in this 
effort are underway. Verification will not be further addressed in this paper. Validation is 


the process of determining the degree to which a model is an accurate representation of 
the real world from the perspective of the intended uses of the model. The last three 
phrases of the validation definition imply that an accuracy requirement has been defined. 
In this paper no accuracy requirement will be defined; our goal is only to determine the 
degree of accuracy. 

Scope 

The scope of this effort is limited to the simulation and evaluation of the results of 
modeling the RCM-1 test case titled “Penn State Preburner Combustor”. The details of 
this test case is described in other documentation and will not be described further here 
unless the explanation of the modeling approach demands otherwise. 

The specific goals of this effort are a) to achieve steady state solutions and 
demonstrate iteration and grid convergence and b) to compare the results of the chamber 
wall heat flux with experimentally determined values and provide a detailed critique of 
the comparison. In the accomplishments of these goals we will consider the effects of 
different turbulence models and of preconditioning and different hybrid grid strategies. 

Computational Tools 

Loci-CHEM is a finite-volume flow solver for generalized grids developed at 
Mississippi State University supported in part by NASA and NSF. Loci-CHEM uses high 
resolution approximate Riemann solvers to solve finite-rate chemically reacting viscous 
turbulent flows. Preconditioning 4 is available for low Mach number applications. 
Various chemical reaction mechanisms are available; the model used in this study for 
hydrogen/oxygen is a 6 species 28 reaction model 5 as shown in Table 1. Thermodynamic 
properties are provided via a standard partition function formulation which calculates the 
specific heats, internal energies and entropies of each individual perfect gas species. 
Several turbulence models are available; the Mentor Shear Stress Transport 6 (SST) and 
original Wilcox 7 k-omega models were used in the current study. Details of the 
numerical formulation are presented in the Loci-CHEM user guide. 8 Loci-CHEM is 
comprised entirely of C and C++ code and is supported on all popular UNIX variants and 
compilers. Parallelism is supplied by the Loci 9 framework which exploits multi-threaded 
and MPI libraries to provide parallel capability. 



Table 1. The six-species hydrogen oxygen reaction model. Units are inks. kf=AT n 

c (-e/rt) ^ determined from Kc. 


Step 

Reaction 

A 

N 

E/R (K) 

Non M-body Reactions 

1 

H20 + 0 <-> 2 OH 

5.8xl0 lu 

0 

9059 

2 

H20 + H <-> OH + H2 

8.4xl0 lu 

0 

10116 

3 

02 + H <-> OH + O 

2.2xlO n 

0 

8455 

4 

H2+0 <-> OH + H 

7.5xl0 lu 

0 

5586 

M-body Read 

tions 

5 

H2 + M <-> 2 H + M 

5.5xlO n 

-1 

51987 

6 

02 + M <-> 2 O + M 

7.2xl0 15 

-1 

59340 

7 

H20 + M <-> OH + H + M 

5.2xl0 18 

-1.5 

59386 

8 

OH + M <-> 0 + H + M 

8.5xl0 15 

-1 

50830 


All grids are generated using GRIDGEN from Pointwise Corporation. 10 Hybrid 
grids were generated by extruding connecters representing the axi-symmetric solid 
surfaces normal to the solid surface into the computational domain with an initial defined 
spacing and stretch rate. The extrusion was terminated at the location where the extruded 
cells possessed an aspect ratio of approximately unity. This allowed a good quality 
transition from structured to unstructured cell types. The remainder of the computational 
domain was gridded using unstructured cells with a boundary decay factor of 
approximately 0.995 to avoid radical coarsening of the unstructured grid. 

Computational Model 

A sketch showing the boundary conditions for all eight cases is shown in Figure 2. 
The inlet walls and post-tip are adiabatic, no-slip walls. The oxygen and hydrogen inlets 
are both modeled as a fixed mass flow rate. Inlet boundary values for turbulence 
quantities, k and omega, were specified as 5x1 O' 6 m 2 /s 2 and 500 1/s respectively. While 
these values are lower than those recommended 6 , the turbulent field develops in the 
injector inlets and previous experience has shown that the injector inlet distances in this 
case allow no sensitivity of the chamber solution to the injector inlet turbulence quantities 
specification. The temperature profile specified in the RCM-1 case definition document is 
imposed as a boundary condition on the no-slip chamber wall. 

The faceplate and nozzle wall temperature boundary conditions are set as isothermal 
fixed temperatures with no-slip conditions for the flow. The faceplate was set to the 
value of the first data point, 639 K, from Table 3 of the RCM-1 case definition document. 
The nozzle wall temperature was set to the value of the last data point, 510 K. 




Chamber Wall Temperature, 
Fixed in Time, 



Figure 2. Sketch of the computational domain indicting location and type of 

boundary conditions. 

The computational domain was initialized to quiescent steam at 780 K and 750 psia to 
start the steady state simulation. The steam provided a high enough initial temperature, 
such that when the reactant streams reached the post-tip region and mixed, self-ignition 
occurred. Loci-CHEM was executed in local time stepping mode with a physical time 
step of lxl O' 4 seconds which was dynamically reduced locally to not exceed a maximum 
CFL number of 1 0,000 or a change in temperature, pressure or density variables of more 
than 10 percent. With these settings, these cases ignited relatively smoothly, i.e., the 
ignition did not cause any sustained reverse flow of propellants into either the fuel or 
oxidizer inlet tubes. 

Solution convergence is evaluated by a combination of residual drop, total mass flow 
convergence in terms of integrated mass at inlet verses integrated mass at outlet, species 
stored mass convergence integrated over the entire domain volume, and temperature and 
pressure iteration history at selected probe point locations throughout the flow-field. In 
addition the chamber wall heat flux was monitored as a function of solution iteration until 
convergence was identified. For the cases presented in this paper three orders of 
magnitude of residual drop are obtained. This is less than observed in other CFD 
programs due to the use of the local time stepping feature documented above in which the 
time step increases as the solution converges, thus tending to increase residuals. 

Computational domain mass conservation was always achieved to with 0.1 percent of 
the total mass flow. Global mass conservation does not necessarily mean that the mass of 
the individual species is not changing with time or iteration. For example, hydrogen and 
oxygen could be forming into water in some part of the chamber without any net change 
in global mass flow rate. Therefore, we track species conservation by integrating the 
stored mass of each individual species within the entire modeled domain and tracking 
that integrated value versus iteration. The species stored masses are also driven down to 
less than 0.1 percent of the total stored mass of that species. Similarly, the probe point 
pressures and temperatures are also driven down to less than 0. 1 percent variation. 

A typical solution requires about 10,000 time steps to achieve steady state 
convergence. For a 500,000 cell grid, using 16 AMD Opteron 246 processors, steady 
state convergence requires approximately 4 days of wall clock time or 1,600 cpu hours. 


As long as the number of cells per processor remains similar to that described above, the 
Loci-CHEM CFD program is nearly perfectly scalable. Computer resources required for 
other grids used in this study can be scaled accordingly. 

A sample flow field result in shown in Figure 3 for the baseline case. The oxidizer has 
mixed and combusted in the first half of the chamber length. The recirculation in the 
forward end of the chamber extends to about one-third of the chamber length. As 
illustrated in Figure 4, the recirculation is composed of a majority of water vapor and a 
minority of hydrogen on a mass basis. 
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Figure 3. Computed temperature distribution, in degrees K, for the baseline case, 


Grid 2. 
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Figure 4. Close-up views of the temperature, oxygen, water vapor and hydrogen 
mass fraction distributions of the baseline case flow field, Grid 2. 


Results 

As mentioned above, the effects of turbulence model selection, preconditioning 
and grid convergence will be assessed for their impact on the predicted chamber wall heat 
flux. The grids used in the study are described first. 

Grid Refinement Studies 

A total of eight different axisymmetric grids were used to assess grid 
convergence, hereafter referred to as Grid 1 through 8. Grid 1 was generated as part of 
an earlier study 11 and consists of 51,030 structured grid cells. The grid spacing along the 
chamber wall is approximately lxl O' 4 inches which corresponds to a first cell y-plus 



value of approximately 1.3. Due to the limitations of structured gridding techniques, the 
first cell y-plus values for other solid surfaces were much greater than unity and resulted 
in unsatisfactory resolution of the turbulent boundary layer. Attempts to resolve all 
boundary layers to a y-plus value of unity resulted in very large cell counts, very large 
cell aspect ratios, and high concentration of cells in locations of little variable gradient. 
Therefore, all subsequent grids were of the hybrid type. Details of Grid 1 are shown in 
Figure 5. 



Figure 5. Details of Grid 1 

Realizing the structured grid limitations, we decided to use a hybrid grid 
approach. In this approach the turbulent boundary layer is resolved on a structured grid 
grown hyperbolically from the solid surface to the estimated outer edge of the turbulent 
boundary layer. The remainder of the computational domain volume is discretized with 
unstructured cells. In this manner the gridding requirements of the turbulence model can 
be satisfied along all solid-surface boundaries. 

The first, baseline hybrid grid was generated with an estimate of appropriate 
resolution of both the turbulent boundary layers and the non-boundary-layer regions of 
the computational domain. Grid 2 achieves first cell y-plus values ranging from 0.1 to 
0.75 for all solid surfaces except the converging-diverging nozzle, where it ranges from 
0.5 to 3.3. The extrusion stretch rate is nominally 1.2 for the structured layer of cells. 
The unstructured potion of the grid is constructed of triangular cells with sides ranging 
from 3xl0' 4 to 1.8xl0' 2 inches. The total number of cells used to construct Grid 2 is 
496,496. Details of Grid 2 are shown in Figure 6. 




Figure 6. Details of Grid 2. 


We addressed the grid refinement of the non-boundary layer region and the 
boundary layer region separately. First, we reduced the maximum triangular cell side 
length from 1 ,8xl0' 2 10 9.0xl0’ 3 inches in the non-boundary layer portion of the grid. The 
resulting grid, Grid 3, was constructed of a total of 743,936 cells. Details of the head end 
of the chamber of Grid 3 are shown in Figure 7. Comparison with Grid 2 in Figure 6 
shows the refinement of the chamber region as a result of the change. 



Figure 7. Details of Grid 3. 

Then we created a sub-region in which the upstream half of the flame is located. 
We specified the maximum triangular cell side length in this sub-region as 6.2x1 O' 3 
inches. The properties for the main unstructured region were unchanged form those used 
in Grid 2. This new grid, Grid 4 was constructed of 782,872 cells. 

We then combined the specifications of the unstructured regions from Grid 3 and 
Grid 4 and created Grid 5, which is constructed of 993,848 cells. A comparison between 


Grid 2 and Grid 5 is contained in Figure 8. The refinement in the flame sub-region 
results in a factor of three better spatial resolutions than Grid 2. 



Figure 8. Comparison of Grid 2 (black) and Grid 5(red) in the near post tip region 

and the downstream region. 


Next, we addressed the grid refinement of the boundary layer region. There are 
three quantities that can be readily changed in a structured extruded grid. They are initial 
spacing normal to the wall, the growth or stretch rate as the cells are extruded from the 
solid boundary and the cell spacing parallel to the solid wall. 

The stretch rate was reduced from 1 .2 in Grid 2 to 1.1 for Grid 6 while holding 
the initial cell spacing nearly constant. Grid 7 increases the initial cell spacing while 
retaining the stretch rate of 1.1. Grid 8 reduces the initial cell spacing both parallel and 
normal to the chamber wall relative to Grid 2 while retaining the stretch rate of 1 . 1 . See 
Figures 9 and 10 for a summary of the main characteristics of the chamber wall boundary 
layer grides used in the grid sensitivity evaluation. 
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Figure 9 Variation of grid initial spacing normal to the chamber wall for all hybrid 

grides. 
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Figure 10. Variation of grid initial spacing parallel to the chamber wall for all 
hybrid grides as a function of perimeter distance from the injector outer diameter 

corner. 


Computational Results 

NASA/MSFC has been studying this RCM-1 case for some time. In 2004, an 
initial study 11 was conducted using Grid 1 with a previous, version of the Loci-CHEM 
CFD program. Although the inlet and boundary conditions are not precisely those of the 
RCM-1 definition, they were a good approximation of the RCM-1 case. The predicted 


heat flux compared to the experimentally determined heat flux is shown in Figure 11. 
The Mentor Baseline and Shear Stress Transport models are a better qualitative and 
quantitative representation of the experimentally determined heat fluxes than the Wilcox 
model. They both indicate the first heat flux peak is the largest magnitude. 

This comparison was repeated using the hybrid structured-unstructured Grid 2 and 
the results are shown in Figure 12. A similar result is obtained with a slight quantitative 
advantage when using the SST model. This perceived advantage is based on the desire 
that a computational prediction to support design must be capable of predicting the 
location and magnitude of the greatest threat environments. The heat flux rise rate is 
important also, but the Wilcox model, which better predicts the rise rate does not do well 
on predicting the peak heat flux details. All models equally over-predict the heat flux 
downstream of the 5 inch location. 



Figure 11. Comparison of heat flux computed with three different turbulence 
models in Loci-CHEM version 2 to the experimentally determined value. 
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Figure 12. Comparison of heat flux computed with three different turbulence 
models to the experimentally determined value. 
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Figure 13. Profiles of normalized chamber near wall values for Grid 2. 


Some discussion of the predicted dual peak nature of the heat flux profile is in 
order. Experimental data does not clearly indicate a dual peak of the heat flux profile. 
To gain some insight into the predicted nature, normalized values of predicted 
temperature, velocity magnitude, turbulent kinetic energy, eddy viscosity and heat flux as 
a function of distance from the faceplate at a radial-location of 0.74 inches, i.e., 0.1 
inches form the chamber wall are shown in Figure 13. The heat flux peaks at about 2.5 
inches downstream of the faceplate. The injector flow attaches to the wall downstream of 
this location at about 3.65 inches downstream of the faceplate as indicated by the velocity 
magnitude disappearing. The velocity magnitude, turbulent kinetic energy and eddy 
viscosity all peak upstream of the heat flux peak. The gas-phase temperature continues to 
rise throughout the first two-thirds of the chamber. The reductions in three of the four 
variables influencing the evolution of the heat flux downstream of about 2 inches cause 
the heat flux to reduce in the downstream direction after its first peak. After the flow 
attachment location, the velocity magnitude outside of the boundary layer increases while 
the turbulent kinetic energy and eddy viscosity continue to decrease slowly. This enables 
a second rise in the heat flux due to increased convective heat transfer. This rise is 
terminated by the plateau of the external boundary gas-phase temperature and the 
continued drop in turbulent kinetic energy and eddy viscosity, which results in a second 
peak in the predicted heat flux profile. 

Next, the effect of preconditioning on the predicted heat flux was assessed using 
Grid 2 and the SST turbulence model. The Loci-CHEM implementation of 
preconditioning requires the specification of a global mach number scale, M®. 
Experience has shown that reasonable accuracy can be obtained for Mach number down 
to about 10 percent of this value. We inspected the flow field computed using Grid 2 and 
consequently set the MMo 0.1, which allows good resolution down to a Mach number of 
about 0.01. The effect of preconditioning on the predicted heat flux is shown in Figure 
14. Although the activation of preconditioning caused no significant change in the 
predicted heat flux, its use did significantly increase the efficiency of obtaining 
convergence. Due to this result, all subsequent computations were executed using 
preconditioning and all figures in this paper using the description baseline are those 
computed using preconditioning. 
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Figure 14. Effect of preconditioning on predicted heat flux for the parameter value 

of Moo of 0.1. 

To asses the effects of grid convergence, predicted heat fluxes using Grids 2 
through 5 are compared in Figure 15. These non-boundary layer region refinements 
systematically decrease the heat flux by approximately five percent in the upstream 
portion of the chamber and increase the heat flux by a similar amount in the downstream 
portion of the chamber. 
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Figure 15. Comparison of predicted heat flux for the series of grides refining non- 
boundary layer portion of the computational domain. 

Refinement of the boundary layer region is represented in Grids 6, 7 and 8. The 
range of the chamber wall first cell y-plus values are shown in Figure 16. Clearly for all 
cases the first cell y-plus values are appropriately small. The effects of the various 
boundary layer refinements on the predicted heat flux are shown in Figure 17. The net 
effect of the boundary layer refinements is to increase the heat flux in the upstream 
portion of the chamber and slightly decrease the heat flux in the downstream portion of 
the chamber. 
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Figurel6. Chamber wall first cell Y+ values for all hybrid grides. 
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Figure 17. Comparison of predicted chamber wall heat fluxes for the series of 
grides used to assess boundary-layer region grid independence. 


In terms of grid refinement, it is clear that the initial grid, Grid 2, is essentially an 
adequate grid for the prediction of chamber wall heat fluxes. The effects of the grid 
refinement in the non boundary-layer verse boundary layer are competing, but small. It is 
the opinion of the authors that in this progression of grides, that grid resolution effects 
have been driven down to the level that they are comparable to the effects of small 
changes in grid quality. 


Conclusions 

Conclusions can be drawn in three specific areas. The first is that in regards to 
the RCM-1 test case, the use or lack of use of preconditioning has no significant effect on 
the accuracy of the predicted chamber wall heat flux. However, the use of 
preconditioning does result in a significant increase in the efficiency of obtaining 
converged solutions to this test case when using the Loci-CHEM CFD program. 

The second area is the sensitivity of the predicted heat transfer to the turbulence 
model chosen. The Menter Shear Stress Transport model provides superior performance, 
relative to the data in the upstream portion of the chamber, for the heat flux rise rate and 
peak value. This point is further amplified by the dissection of the computed solution 
which illustrates how the predicted twin peak nature of the heat flux arises from the 
distributions of the turbulent quantities relative to the velocity and thermal fields. 

The last area is that reasonable grid independence has been demonstrated. The 
initial grid, Grid 2 with about 500,000 cells has been proven to be essentially providing 
the grid independent answer. 

Overall, this effort demonstrates the limit of the predictive ability of the Loci- 
CHEM CFD program on the RCM-1 test case. Predictions of the initial heat flux rise rate 
and heat flux peak are satisfactory for engineering use. The heat flux is significantly 
over-predicted in the downstream portion of the chamber, rendering these predictions 
questionable for engineering use in this region. 

Future Work 

The mixed results found in this work illustrate the need to significantly improve 
the predictive ability of CFD methods. The demonstrated nature of the balances that lead 
the specific behavior of the heat flux prediction indicate that turbulence modeling should 
be addressed to achieve this improvement. 

Improvements can be sought by investigating and improving RANS turbulence 
models, investing the use of LES methods and/or the use of hybrid RANS-LES methods. 
In any case, the evaluations should follow a very careful assessment of the proposed 
improvement via a hierarchy 12 of problems ranging from extremely simple and 
systematically progressing to more complex problems. 

The grid refinement results presented in this paper can be further evaluated by 
applying advanced techniques 13 to the assessment of solution convergence. We will be 
pursuing these improvements in the future. 
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